function result=canny(image, T1, T2, Sigma)
    image = im2double(image);

    if ~exist('T1','var')
        T1 = 0.7;
    end
    if ~exist('T2','var')
        T2 = 0.1;
    end
    if ~exist('Sigma','var')
        Sigma = 0.5;
    end

    filter=fspecial('gaussian', [3 3], Sigma);
    r = imfilter(image, filter);
    
    rx = edge(r, 'sobel', 0.25, 'horizontal');
    ry = edge(r, 'sobel', 0.25, 'vertical');
    
    rx = rx .* rx;
    ry = ry .* ry;
    r = sqrt(rx + ry);
    
    s = size(image);
    angles = zeros(s(1), s(2));
    suppression = zeros(s(1), s(2));

    for y = 2:s(1)-1
        for x = 2:s(2)-1
            % ------------------------ Define direction to follow
            
            a = atand(ry(y, x) / rx(y, x));
            if a < 0
                a = a + 360;
            end
            if a >= 360
                a = a - 360;
            end
            if a >= 180
                a = a - 180;
            end
            if a < 22.5
                a = 0;
            else if a < 67.5
                    a = 45;
                else if a < 112.5
                        a = 90;
                    else if a < 157.5
                            a = 135;
                        else
                            a = -1;
                        end
                    end
                end
            end
            angles(y, x) = a;
            
            % ------------------------ Non maximum suppression
            
            if a == 0
                if (r(y, x) >= r(y, x-1) && r(y, x) >= r(y, x+1))
                    suppression(y,x) = r(y,x);
                else
                    suppression(y,x) = 0;
                end
            else if a == 45
                    if (r(y, x) >= r(y+1, x-1) && r(y, x) >= r(y-1, x+1))
                        suppression(y,x) = r(y,x);
                    else
                        suppression(y,x) = 0;
                    end
                else if a == 90
                        if (r(y, x) >= r(y-1, x) && r(y, x) >= r(y+1, x))
                            suppression(y,x) = r(y,x);
                        else
                            suppression(y,x) = 0;
                        end
                    else if a == 135
                            if (r(y, x) >= r(y-1, x-1) && r(y, x) >= r(y+1, x+1))
                                suppression(y,x) = r(y,x);
                            else
                                suppression(y,x) = 0;
                            end
                        end
                    end
                end
            end
        end
    end

    % ------------------------ Hysteresis treshold
    
    result = zeros(s(1), s(2));
    changes = 1;
    
    for x = 2 : 2 : s(2)-1
        for y = 2 : s(1)-1
            T1sum = 0;
            if result(y-1,x-1) > T1
                T1sum = T1sum + 1;
            end
            if result(y,x-1) > T1
                T1sum = T1sum + 1;
            end
            if result(y+1,x-1) > T1
                T1sum = T1sum + 1;
            end
            if result(y+1,x) > T1
                T1sum = T1sum + 1;
            end
            if result(y+1,x+1) > T1
                T1sum = T1sum + 1;
            end
            if result(y,x+1) > T1
                T1sum = T1sum + 1;
            end
            if result(y-1,x+1) > T1
                T1sum = T1sum + 1;
            end
            if result(y-1,x) > T1
                T1sum = T1sum + 1;
            end
            if suppression(y, x) > T1 && T1sum <= 3
                result(y, x) = suppression(y, x);
            end
        end
    end
    
    for x = 3 : 2 : s(2)-1
        for yy = 2 : s(1)-1
            y = s(1) - yy + 1;
            T1sum = 0;
            if result(y-1,x-1) > T1
                T1sum = T1sum + 1;
            end
            if result(y,x-1) > T1
                T1sum = T1sum + 1;
            end
            if result(y+1,x-1) > T1
                T1sum = T1sum + 1;
            end
            if result(y+1,x) > T1
                T1sum = T1sum + 1;
            end
            if result(y+1,x+1) > T1
                T1sum = T1sum + 1;
            end
            if result(y,x+1) > T1
                T1sum = T1sum + 1;
            end
            if result(y-1,x+1) > T1
                T1sum = T1sum + 1;
            end
            if result(y-1,x) > T1
                T1sum = T1sum + 1;
            end
            if suppression(y, x) > T1 && T1sum <= 3
                result(y, x) = suppression(y, x);
            end
        end
    end
    
    while changes ~= 0
        changes = 0;
        for x = 2 : s(2)-1
            for y = 2 : s(1)-1
                sum = 0;
                if result(y-1,x-1) > 0 
                    sum = sum + 1;
                end
                if result(y,x-1) > 0 
                    sum = sum + 1;
                end
                if result(y+1,x-1) > 0 
                    sum = sum + 1;
                end
                if result(y+1,x) > 0 
                    sum = sum + 1;
                end
                if result(y+1,x+1) > 0 
                    sum = sum + 1;
                end
                if result(y,x+1) > 0 
                    sum = sum + 1;
                end
                if result(y-1,x+1) > 0 
                    sum = sum + 1;
                end
                if result(y-1,x) > 0 
                    sum = sum + 1;
                end
                if result(y, x) == 0 && suppression(y, x) > T2 && sum <= 3 && sum > 0
                    result(y, x) = suppression(y, x);
                    changes = changes + 1;
                end
            end
        end
    end
end